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Urchin: A Reverse Ray Tracer for Astrophysical Applications 
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ABSTRACT 

We describe URCHIN, a reverse ray tracing radiative transfer scheme optimised to model 
self-shielding from the post-reionisation ultraviolet (UV) background in cosmological sim- 
ulations. The reverse ray tracing strategy provides several benefits over forward ray tracing 
codes including: (1) the preservation of adaptive density field resolution (2) completely uni- 
form sampling of gas elements by rays; (3) the preservation of galilean invariance; (4) the 
ability to sample the UV background spectrum with hundreds of frequency bins; and (5) ex- 
act preservation of the input UV background spectrum and amplitude in optically thin gas. 
The implementation described here focuses on Smoothed Particle Hydrodynamics (SPH). 
However, the method can be applied to any density field representation in which resolution 
elements admit ray intersection tests and can be associated with optical depths. We charac- 
terise the errors in our implementation in stages beginning with comparison to known analytic 
solutions and ending with a realistic model of the z = 3 cosmological UV background inci- 
dent onto a suite of spherically symmetric models of gaseous galactic halos. 
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1 INTRODUCTION 



Cosmological gas dynamics simulations (e.g. |Cen et al.|1994}|The-| 
|uns et af||1998b|a[ Hernquist et al. 1996), semi-analytic models 
(e.g. Bi & Davidsen 1997) and analytic calculations (e.g. Schaye 
|2001[ l have enabled us to understand the connection between galax- 
ies, the intergalactic medium (IGM), and the large-scale structure 
of the Universe. HI Lyman-a absorbers in the spectra of distant 
quasars are a particularly useful observational probe of this struc- 
ture. Models indicate that at redshifts z = 2 — 4, the majority 
of Lyman-a forest absorption lines (i.e. lines with column densi- 



ties TVhi < 



) arise in filamentary structures in the 



IGM with the highest column density lines being associated with 
the circumgalactic medium of galaxies. These lines are produced 
in systems which are highly ionized by the ultraviolet (UV) back- 
ground (see Meiksin 2009 . for a recent review). An absorber with 
a column density iVm > 10 17 ' 2 cm -2 has an optical depth greater 
than unity for Lyman-Limit photons and is called a Lyman-Limit 
System (LLS). Above a column density of JV H i = 10 20,3 cm~ 2 , 
'damping wings' due to the natural line-broadening of the Lyman- 
a line are detectable, and the system is called a damped Lyman-a 
absorber (DLA). These absorbers probe the interface between the 
IGM and galaxies as well as the interstellar medium (ISM) of the 
galaxies themselves. Hydrogen begins to self-shield from the UV 
background in the LLS column density range, and the reduction 
in ionising flux plays a major role in setting the ionisation state of 
these absorbers. The URCHIN code described in this paper is de- 
signed to model HI self-shielding in the LLS and DLA range. 



At present, the largest observational catalogues of self- 
shielded absorbers are produced through semi-automated searches 
(e.g. |Prochaska & Herbert-Fort|2004"{ |Noterdaeme et al.|2012| > of 
data from the Sloan Digital Sky Survey (SDSff}. Current and 
planned expansions to the SDSS such as the Baryon Oscillation 
Spectroscopic Survey (BOSS, [Schlegel et al.|2007l ) and BigBOSS 
(Myers et a l.|2012) will increase the amount of available data by 
a factor of ten. Due to atmospheric absorption of the rest frame 
Lyman-a transition, ground based surveys for DLAs are limited 
to redshifts z > 1.6. Surveys for LLSs require spectral cover- 
age of the Lyman-Limit transition for an accurate determination 
of Vhi and are therefore limited to redshifts of z > 2.5 when per- 
formed from the ground. The new Cosmic Origins Spectrograph?] 
on the Hubble Space Telescope provides significant capacity to 
probe lower redshift systems (e.g. Battisti et al. 20121 while the 
Advanced Camera for Surveys and Wide Field Camera 3 have re- 
cently been used to complete a survey for LLSs in the redshift range 
1.0 < 2 < 2.6 ( [O'Meara et al.|2012| >. 

In addition to Lyman series transitions, post-reionisation neu- 
tral hydrogen can also be effectively probed using the 21 -cm emis- 
sion line. The most recent determination of the local HI mass func- 
tion is from the Ar ecibo Legacy Fast ALFA (ALFALFA) survey 
I Martin et al.|2010| which will have detected r;3x 10 4 galaxies 
in HI 21-cm out to z = 0.06 when it is complete. The Square Kilo- 
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meter Array (SKjArb represents the long term future for this type 
of radio astronomy however construction will not begin for sev- 
eral years. In preparation of this, a host of pathfinding telescopes 
( ASKAlQ MeerKAlQ WSRlQ EVLA[] ) will soon make 21- 
cm emission surveys as data rich as their optical and near-infrared 
counterparts. In addition, pilot surveys for 21 -cm absorption in the 
spectra of radio bright sources have shown potential (e.g. Gupta 
|et al.|20 10; Darling et al. 201 1 ) and at least two of these pathfinders 
(ASKAP and MeerKAT) will also perform large, blind, absorption 
surveys. 

The combined output from these surveys will provide transfor- 
mative information on the gas content of galaxies and their modes 
of accretion. It will also generate samples that trace the large scale 
structure of the Universe with different biases than those of op- 
tically selected samples. A prerequisite for making model predic- 
tions of hydrogen emission or absorption is the accurate calculation 
of the distribution of neutral hydrogen (e.g. Duffy et al. 2012; van 
|de Voort et al.|2012) . Methods that accomplish this task with a min- 
imum of free parameters will be able to take full advantage of ob- 
servational data. These issues motivate many of the design choices 
for Urchin . 

The standard approach of treating the post-reionization UV 
background in cosmological simulations is to impose a spatially 
uniform but time varying radiation field, and calculate the Hi frac- 
tion in the optically thin limit. Pioneering work on modeling self- 
shielding in gas dynamics simulations was done by Katz et ah] 
( 1996) and Haehnelt et al. ( 1998) by post-processing column den- 
sity maps. More recent theoretical work has incorporated radiative 
transfer through 3-D density fields to calculate the attenuation of 
the UV background in dense gas (Razoumov et al. 2006; Kohler & 
Gnedin 2007, Po ntzen et al.|2008||Altay et al.|201 l||Mc Quinn et al. 
2011[|Yajima et al.|2011||Fumagalli et al.|2011||Cen|2012||Erkal| 



et al.|2012| >. 

In this paper we present and test URCHIN, a reverse ray trac- 
ing scheme designed to calculate self-shielding corrections in the 
post-reionisation Universe. The code can be applied to any den- 
sity field representation (e.g. particles, adaptive grids, unstructured 
meshes) in which the resolution elements can be associated with 
optical depths and subjected to ray intersection tests. The main 
benefits of URCHIN are: (1) preservation of the adaptive density 
field resolution present in many gas dynamics codes; (2) uniform 
sampling of gas resolution elements with rays; (3) preservation of 
galilean invariance; (4) high spectral resolution; and (5) preserva- 
tion of the standard uniform UV background in optically thin gas. 
The format of this paper is as follows. In §2 we introduce our nota- 
tion and review some basic physics related to radiative transfer. In 
§3 we give a general description of our reverse ray tracing approach 
and place it in context by comparing it to alternative approaches. In 
§4 we discuss the details of our implementation using smoothed 
particle hydrodynamics (SPH) density fields. In §5 we present the 
results of tests meant to validate URCHIN and in §6 we discuss the 
results, suggest improvements for future versions of the code, and 
conclude. 
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2 DEFINITIONS AND BASIC PHYSICS 

In this section, we define our notation and review some of the rele- 
vant physics. All simulations discussed in this work utilize a cubic 
simulation volume. For brevity, we will refer to any of these sim- 
ulation volumes as boxes and their six faces as walls. When refer- 
ing to distances, we will distinguish between proper and comoving 
measures using the prefixes 'p 'and 'c '(e.g. pkpc, cMpc). In this 
work, we consider only hydrogen and leave the inclusion of other 
elements, particularly helium, to future work. For our description 
of the radiation field, we adopt the notation of Rybicki & Lightman 
(1986). 

The specific intensity, /„ = dE/dA dfi dt dv, fully charac- 
terises a radiation field and is defined as the energy dE passing 
through an area element dA into a solid angle element dQ in time 
dt due to photons with frequency between v and v + dv. Several 
useful characterisations of the radiation field can be expressed as 
integrals over this quantity. Considering photons with frequency 
fth < v < qvth, and an optically thin medium, we can write the 
number density of hydrogen ionising photons and the photoionisa- 
tion rate at the point x respectively as: 



n 7 (x) 

r(x) 



h 



dfdil, 



hv 



hv 
dv dQ. 
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where hv t h = E t h and a are the ionisation energy and pho- 
toionisation cross-section of hydrogen. In the case of a medium 
with finite opacity, the above quantities can be written in terms of 
the optically thin value of I v by making the replacement I v — > 
/„ exp [— r(v, fl)] where r is the optical depth between the sources 
producing I v and x. 

The frequency averaged ('grey') photoionisation cross section 
is defined as o" grcy = T (cn^) -1 . Under the grey approximation, 
every polychromatic spectrum (characterised by n 7 and F) corre- 
sponds to an equivalent monochromatic spectrum with the same 
n 7 and T. The flux of the monochromatic spectrum is fixed by n 7 
and the energy of the photons in the monochromatic spectrum is 



E u 



hv gvcy with v gvcy implicitly defined by a(i/ gIcy ) = a g 



3 REVERSE RAY TRACING METHOD 

URCHIN is designed to efficiently model the residual neutral hydro- 
gen in the post-reionisation universe. In this section, we describe 
the algorithms used in URCHIN and conceptual departures from al- 
ternative radiative transfer codes. We begin with a brief description 
of the standard treatment of the UV background in cosmological 
gas dynamics simulations. 



3.1 Standard Treatment of Post-Reionisation UV 
Background 

In the post-reionisation Universe, cosmic hydrogen is kept highly 
ionised by a pervasive UV background (Gunn & Peterson 1965) 
produced by galaxies and quasars. Quantitatively, the volume av- 
eraged neutral frac tion x = n HT /n H for redshifts z ^ 6 is on 
the order of 10~ 4 (Fan et al.|2006| JBecker et aL]|2007). A rapid 
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transition to higher neutral fractions signals the end of reionisation, 
evidence for which has recently been observed in the form of a 
damping wing in the spectrum of a z ~ 7 quasar (Mortlock et al. 
|2011[>. At lower redshifts, the UV background determines both the 
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ionisation state and temperature of gas in the IGM and sets the rate 
at which denser gas can cool, accrete onto small galaxies, and form 
stars jEfstathiou|1992l|Okamoto et al.|2008| . 

The most widely used treatment of the post-reionisation UV 
background in cosmological simulations, is based on three approx- 
imations: (1) optically thin gas; (2) a spatially uniform UV back- 
ground; and (3) photo/collisional ionisation equilibrium. These ap- 
proximations are valid for the majority of cosmic gas, i.e. for highly 
ionised hydrogen; however, they break down for the majority of 
gas observable in HI surveys. The approximation of optically thin 
gas begins to break down in absorption systems that probe accre- 
tion (LLSs) from the IGM onto galaxies and completely fails in 
regions of significant self-shielding (DLAs) where the strongest HI 
signals arise. This approximation is the most important to remedy 
for HI surveys. The second approximation involves disregarding 
large scale gradients in the UV background as well as point sources. 
These fluctuations can have an effect on absorber statistics (e.g. 
|Rahmati et al.|2013[[Wjima et al.|20TT]|Croft|2004| > but the mag- 
nitude is not as large as that due to self-shielding. A notable excep- 
tion is absorbers in the proximity zones of bright sources (Schaye 
2006). The third approximation involving photo/collisional equi- 
librium likely holds for dense gas where HI is self-shielded and 
recombination times are short compared to UV background vari- 
ability, but will break down near variable sources and in less dense 
gas. Relaxing these approximations is an efficient way to improve 
model predictions for HI surveys. In this release of URCHIN we 
focus on self-shielding. 



3.2 Reverse Ray Tracing - Motivation 

Numerical techniques for continuum radiative transfer have been 
developed based around ray tracing methods, (e.g. Nakamoto et al. 
200T] [Maselli et al.|2003||Razoumov & Cardall|2005HSusa|2006| 
Whalen & Norm an 2006, Altay et al. 2008), the closely related 
method of characteristics (e.g. Mellema et al. 2006; Rijkhorst et al. 
2006), angular moments of the radiative transfer equation (e.g. 
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2009, Finlator et al. 2009), and transport on unstructured meshes 
(e.g. [pawlik & Schaye||2008l |20T7| |Paardekooper et al.||20T0) . 
A detailed comparison between many codes currently in use is 
documented in the Cosmological Radiative Transfer Comparison 
Project piiev et al.|2006l[2009| >. 

The above methods arc all based on following radiation from 
its source to the point where it is absorbed or scattered. When deal- 
ing with the post-reionisation UV background, the goal is to build 
up a radiation field that is known from both theoretical and ob- 
servation studies to be mostly uniform. In these methods, the UV 
background field at a given point is the sum of radiation that has 
been transported from all sources being considered. However, in 
many methods this is true only in a statistical sense. For exam- 
ple, in Monte Carlo ray tracers such as SPHRAY, resolution ele- 
ments are updated whenever a ray intersects them. The ray only 
carries information about the flux from one source, but if enough 
rays are traced, each resolution element will be updated by the rays 
from many sources. In URCHIN we attempt a different type of so- 
lution to this problem. As opposed to building up a mostly uni- 
form UV background from multiple sources, we begin with the 
standard approximations described above (optically thin gas, uni- 
form field, photo-collisional equilibrium) and then calculate devia- 
tions from it. For the majority of cosmic gas, the optically thin pho- 
toionisation rate from this uniform field, r thm , is in fact a good ap- 
proximation. The current version of URCHIN relaxes the optically 



thin approximation from the standard treatment by attenuating this 
uniform radiation field in denser regions that begin to self-shield. 
The fraction of gas (by mass or volume) where this correction is 
necessary is guaranteed to be small by the nature of the problem 
allowing us to concentrate the available computational resources 
where they are needed. In future versions of the code we will re- 
lax the second and third approximations in the standard treatment 
by adding large scale gradients, proximity regions, and consider- 
ing non-equilibrium effects. Currently URCHIN operates on static 
density fields as a post-processing step, but in principle could be 
coupled to existing gas dynamics codes in a straight-forward way 
(e.g. in a framework such as described in Portegies Zwar t et al.| 
[2009l >. 



3.3 Reverse Ray Tracing - Algorithm 

Consider a density field discretised into resolution elements labeled 
i — 1, • • • , N. We will refer to these resolution elements as parti- 
cles; however, our reverse ray tracing technique can be applied to 
any density field representation in which the resolution elements 
can be associated with optical depths and subjected to ray inter- 
section tests. From each particle, we cast -/V ray ray segments out 
to a distance Z ray , in directions that uniformly cover the solid an- 
gle around the particle. We use the HEALPix algorithm (Gorski 
|et al.|20 05) to determine ray directions. We then calculate the op- 
tical depth Tfc along each of these ray segments and sum over all 
rays to obtain an estimate of the self-shielded photoionisation rate, 



^ r ln , at the location of each particle, 
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(3) 



For most particles the optically thin approximation is very good and 
the effective optical depth will be small, r e g <S 1, along all N lay 
directions. However, for the small fraction of particles that dom- 
inate the HI abundance, the effective optical depth will be large 
(i"cff 3> 1). We then use r ahld in an analytic solution for the equi- 
librium neutral fraction x (Eq. |A9| in the Appendix) to update the 
ionisation state of the particle. An altered neutral fraction for one 
particle leads to altered optical depths along ray segments that pass 
through it, and hence we iterate this procedure until the neutral frac- 
tion converges for all particles. This typically happens in tens of it- 
erations with those few particles at the threshold between optically 
thin and optically thick converging last. l Tay and N Tay are numerical 
parameters of the scheme, and we have found that in fully cosmo- 
logical runs at z ~ 3, values of i ray = 100 pkpc and -/V ray = 12 
lead to converged column density distribution functions. 



3.4 Reverse Ray Tracing - Advantages 

Reverse ray tracing presents some important advantages in the post- 
reionization regime, related to the dramatic differences in the char- 
acter of the radiation field during and after reionization. In the fol- 
lowing sub-sections, we focus on questions that typically occur in 
ray tracing schemes attempting to model the post-reionization cos- 
mological UV background: 1) Where should rays originate and ter- 
minate? 2) How should rays be traced to uniformly sample the gas 
elements? 3) How much spectral resolution can be attained? and 4) 
What optimizations can be applied? We will show that answers to 
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these questions are simpler in reverse ray tracing schemes than in 
forward versions. 



3. 4. 1 Where Should Rays Originate ? 

In forward ray tracing schemes, the UV background at a given point 
in space is the result of transporting photons along rays originating 
at sources of radiation. These sources can be divided into two cate- 
gories: those inside the box (internal sources) and those outside the 
box (external sources). Modelling the UV background from first 
principles using only internal sources imposes stringent conditions 
on the size of the box. A reasonable absolute minimum scale is one 
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mean free path at the Lyman limit, AJ^|\ 
find that this mean free path depends on redshift z approximately 
as X(z) « [(48 ± 2.1) - (38.0 ± 5.3) (2 - 3.6)] h^ pMpc for 
3.6 < z < 4.3. At z — 3.6 this is larger than the majority of 
cosmological gas dynamics simulations able to resolve any galaxy 
formation processes and AJ^^ becomes larger at lower redshifts. 
For this reason, forward ray tracing schemes inevitably must rely 
on tracing many rays from external sources in addition to internal 
ones. 

The contribution from external sources is usually modeled by 
casting rays inward from the walls. These rays are then followed 
until a given fraction of their initial photon content is absorbed. The 
most straight forward method of choosing ray origins is a pseudo- 
random sampling of points on the walls. The appropriate flux is 
usually determined by requiring that the photoionisation rate in op- 
tically thin gas match an input value. This method has two draw- 
backs. The first, discussed in the next sub-section, is the difficulty 
in uniformly sampling an adaptive density field. The second is the 
artificial gradient in photoionisation rate that will develop between 
regions near the centre of the box, where the background rays will 
be most attenuated, and regions near the walls, where the back- 
ground will be at its optically thin strength. This is what we term the 
galilean invariance problem. This adds ambiguity to any calibration 
of flux from the walls and makes the calibration dependent on the 
box size. In general, tracing rays from the walls leads unavoidably 
to a loss of symmetry in the UV background and difficulty in flux 
calibration. These problems can be alleviated somewhat by tracing 
rays from randomly selected under-dense regions within the box 
(Maselli & Ferrara 2005 ), but not eliminated entirely. Some authors 
(e.g. McQuinn et al. 20 11 1 have circumvented this problem by ex- 



cising small regions around halos and ray tracing them individually 
however this neglects filamentary gas between halos. 

In a reverse ray tracing method, these problems are not 
present. All symmetries of the UV background are maintained and 
the photoionisation rate at a given particle position is independent 
of any properties of the box. In addition, the intensity of the UV 
background in optically thin gas is equal to its value in the standard 
approximation (optically thin gas, uniform field, photo-collisional 
equilibrium) by construction. This conveniently allows for direct 
comparisons between results with and without radiative transfer 
with no need for calibrating fluxes. 



3.4.2 The Uniform Sampling Problem 

Modern gas dynamical simulations discretise density fields with 
adaptive resolution elements. This typically leads to better spatial 
resolution (i.e. smaller resolution elements) in regions of higher 
gas density. A radiative transfer approach which relies on rays uni- 
formly cast from the simulation walls leads to poor sampling of the 



resolution elements. Specifically, the overdense self-shielded re- 
gions where radiative transfer is most important are undersampled, 
and the underdense optically thin regions where radiative transfer 
is unnecessary are oversampled. A scheme in which rays split and 
merge (Wise & Abel 201 1 ) can maintain a constant number of rays 
intersecting each resolution element but comes at the cost of in- 
creased overhead. On the other hand, the reverse ray-tracing algo- 
rithm samples each particle with the same number of rays by con- 
struction. This helps focus computational power where it is needed. 

3.4.3 Spectral Resolution 

In a forward ray tracing scheme, the radiation field at a particle 
position is the summation of contributions from different rays. For 
monochromatic spectra, all attenuation information can be charac- 
terised by the number of photons in a ray. To accurately handle 
multi-frequency spectra, one must substantially increase the num- 
ber of monochromatic rays traced or carry spectral information 
along with each ray. Similarly for moment methods, each frequency 
group must be treated separately by the solver. Low fidelity sam- 
pling of the UV background spectrum can lead to errors up to an 
order of magnitude in HII region, see e.g. (Mirochaet al. 2012) and 
Figure[3] below. 

In reverse ray tracing, the rays simply sample the optical depth 
along a particular direction. Subsequently, arbitrarily complex in- 
put spectra can be numerically integrated over frequency to deter- 
mine a shielded photoionisation rate, and update the neutral frac- 
tion of a particle based on full knowledge of the amplitude and 
spectral shape of the local radiation field. When using 100 fre- 
quency bins in URCHIN, the numerical integration over the UV 
background spectrum consumes a negligible fraction of the com- 
puting time. This allows for studies of spectral hardness and eases 
the inclusion of sharp features such as recombination lines into the 
UV background spectrum. As an added bonus, knowledge of the 
optical depth in all directions around a particle during the update 
allows one to choose the appropriate 'case A' (recombinations to 
all levels) or 'case B' (recombinations to all but the ground state) 
recombination rates in the on-the-spot approximation. 

3.4.4 Optimisations 

The reverse ray tracing approach also lends itself to some important 
optimisations that are not easy to implement in forward ray tracing 
techniques. The most expedient is simply avoiding radiative trans- 
fer where it is unnecessary. As discussed above, particles that are 
not self-shielded and not in proximity zones are treated correctly 
in the uniform and optically thin limit. When looping over the par- 
ticles in reverse ray tracing schemes, those that satisfy this crite- 
rion can be skipped. Identifying these particles is implementation 
dependent; however, in cosmological simulations, the majority of 
post-reionization gas does not need to participate in radiative trans- 
fer. We will discuss our technique for identifying these particles in 
the next section. 

There is also a sense of ray locality and ray independence in- 
herent in the method that can be useful in optimisations. Any sub- 
volume in a box can be treated independent of any other as long 
as a buffer of Z ray is included around the sub- volumes. In addition, 
each ray can be made independent during a single iteration by wait- 
ing to update the neutral fractions of particles until shielded pho- 
toionisation rates have been calculated for all of them. This strategy 
may increase the number of iterations necessary for convergence 
but would make each iteration much faster. 
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4 IMPLEMENTATION DETAILS 

URCHIN applies a specified radiation field to a given density field to 
determine the properties of self-shielded regions. Spatially uniform 
models of the cosmological UV background whi ch provide I v (z), 
are publically available (for example Haardt & Madau 2012 r The 



main parameters in URCHIN are the frequency range of ionising 
photons included in the calculation, the number of rays casts per 
particle, iV ray , their proper length, l la , y , the criterion for choosing 
between type A and type B recombination rates, and the condi- 
tion for convergence 5x to \. The choice of spectrum and frequency 
range determines the optically thin photoionisation rate, p*" 111 , 
from Eq. J2l. Our default runs use u t h < v < 4fth, N Ta . y — 12 
and Z ray = 100 pkpc, which give converged numerical results at 
redshift z = 3. In addition, we use case B recombination rates for 
particles with t«b ^ 1. 

Initially the neutral fraction x of each particle is set to its opti- 
cally thin value, Xthin = njn/flH. This allows for the calculation 
of the HI column density along each ray cast from a particle. Next 
we calculate a shielded photoionisation rate, r shld ; for each parti- 
cle and an effective optical depth, r 8 fl using Eq. p). Particles with 
Toff < T"cft P maintain x = Xthin for all subsequent iterations while 
all other particles are updated each iteration. We continue to loop 
over particles which initially had T e g > r^ s lp until convergence in 
the neutral fraction, |<$a;|/a; < <5xtoi- Typical values for these pa- 
rameters are T c s ff lp = 10 -4 , Sx to i = 10 -3 , l Ta , y = 100 pkpc and 
N ra .y = 12. For these choices we find that more than 99 percent 
of particles have converged neutral fractions after five iterations in 
a cosmological run, with the remaining 1 percent requiring tens 
of iterations. In our current implementation, we update the neutral 
fraction of each particle as soon as its self-shielded photoionisation 
rate has been computed. This update scheme works well provided 
the list of particles is traversed from high to low density, however, 
there is no fundamental restriction on when the particles should be 
updated. Implementations which calculate r shld for each particle 
before performing any updates of the neutral fractions are indepen- 
dent of the order in which the particles are looped over, which could 
be useful in parallel strategies. 



4.1 Reverse Ray Tracing in SPH 

Although URCHIN can be applied to many types of density field 
discretisations (particles, grids, unstructured meshes), in this sec- 
tion we discuss our implementation in Smoothed Particle Hydro- 
dynamics (SPH, |Gingold & Monaghan|1977[[Lucy|1977^ . 



4.1.1 Column Densities in SPH 

Our calculation of SPH column densities makes use of several im- 
provements over the algorithm described in Altay et al. (2008). In 
the SPH formalism the number density of HI at any location r can 
be calculated using the scatter approach as 



ram(r) 



£ 



W(<& 



(4) 



where m p is the mass of a proton, rrn is the mass of particle i in 
hydrogen, Xi is its neutral fraction, qi = |r — Ti\/hi is the dis- 
tance between the particle and and the point r in units of the par- 
ticle's smoothing length hi, and W is the SPH smoothing kernel. 
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The column density through an SPH distribution along a path r(Z) 
parameterised by < I < L can then be written as: 
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(5) 



where the summation is over particles with smoothing volumes in- 
tersected by the path r(Z) and the subscripts in the variable qu in- 
dicate that it is a function of the summation index i and distance 
along the path I. In this way, the calculation of optical depths is 
reduced to calculating which particles are intersected by a ray and 
line integrals through the smoothing kernel W. 

The Gaussian function has many properties that make it a nat- 
ural choice for the smoothing kernel, however, its lack of compact 
support leads to an impractical sum over all particles in the simula- 
tion volume. To remedy this, many SPH codes make use of spline 
functions with an approximately Gaussian shape, such as the Mi 
cubic spline (Monaghan & Lattanzio 1985), 



M t (q) 




6q 2 + 6q i for < q < \ 
for | < q < 1 

otherwise . 



(6) 



In URCHIN we use a truncated Gaussian kernel that allows us to 
obtain line integrals at a given impact parameter analytically. A 
normalised Gaussian function centred at the origin is given by 
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where A 2 = (2cr 2 ) -1 and a 2 is the variance. We truncate G at 
r — h and further demand that G(0, a 2 ) — M±(0, h) to obtain the 
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and the error function is defined as, 



erf(i) = -= / exp(-t 2 )dt. 
\Ar Jo 



(9) 



(10) 



The column density through such a Gaussian kernel at impact pa- 
rameter b, between the limits z\ and zi is then 



I(h,b, zi,Zn) 



/Vexp(-A 2 6 2 



[erf(Az 2 ) - erf(Azi)] . (11) 



We use this kernel to calculate the HI column densities along each 
ray using Eq. pi. Incidentally, the fact that the Gaussian kernel 
can be decomposed into three 1-D functions makes it useful for 
smoothing SPH particles onto 2-D or 3-D grids. 
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4.1.2 Self-contribution to self-shielding in SPH 

When passing from a continuous density field representation to a 
discrete one, Eq. ([A9j for the equilibrium neutral fraction at the 
point r, x — x(r,T,n H ,T,y) becomes partially implicit: in the 
discrete case, F has a dependence on x through the particle's own 
contribution to the optical depth along a ray. This is a generic prob- 
lem of discretisation and has been discussed previously in Abel 
|eTaLl ( fT999) and |Mellema"eTaI1 ( [2r306l l. 

However when the density field is represented with SPH par- 
ticles, further complications arise due to the weighted sum over 
neighbours. For example the HI density at location r depends on the 
neutral fractions of all particles that appear in the sum of Eq. 0), 
This makes a straightforward scheme that updates neutral fractions 
based on r shld calculated at the center of an SPH particle unsta- 
ble. The reason for this instability is illustrated in the top part of 
FigurefT] Consider three SPH particles which are initially optically 
thin but are located in a region that will eventually become self- 
shield. Rays traced outward (in any direction) from each particle 
will probe optical depth contributed by all three. This will in turn 
cause the neutral fraction of each particle to be increased. In the 
next iteration, each ray will find increased optical depth and each 
neutral fraction will increase again. This process will continue and 
eventually spread to adjacent particles causing unphysical growth 
of neutral regions. The cause of this instability is the multi- valued 
relationship between points in space and resolution elements in 
SPH. As a counterexample, consider the uniform grid density field 
represented in the bottom part of Figure [T] This discretisation has 
a single-valued relationship between points in space and resolution 
elements. In this case, changes in the neutral fraction of each el- 
ement do not necessarily change the optical depth encountered by 
each ray. For example, an increase in the neutral fraction of element 
1 has no effect on the optical depth encountered by rays traced from 
elements 2 or 3. This prevents the instability from occuring. 

We resolved this numerical artefact in SPH as follows. Parti- 
cles intersected by each ray are split into two distinct groups labeled 
near and far, with r = T ncar + r far . Near particles of particle i are 
those that contribute to i's (neutral) density in Eq. (El, i.e. the par- 
ticle's neighbours. All other particles that contribute to the sum in 
Eq. J5} are labelled 'far'. The numerical instability only involves 
near particles, therefore we can treat the far particles with the al- 
gorithm described in §3. We model the near particles as a uniform 
slab with the same temperature T and density nu as particle i. The 
thickness of the slab is quantified by calculating the total hydrogen 
optical depth of the near particles at the Lyman Limit defined as 
T near = ATg^Vth- The column density A^ car is calculated as in 
Eq. (pi except all of the neutral fractions are set to unity. 

The radiation incident onto the slab is determined by the user 
supplied spectrum and the optical depth through the far particles, 
/„ exp(— r ar ). For each ray k — 1, • • • , N Tas this incoming flux 
provides a photoionisation rate at the surface of the slab of F fc ar . 
We then solve for the ionisation structure in the slab and associate 
the photoionisation rate at the bottom, r| hld , with the contribution 
from ray k to the total shielded photoionisation rate r shld of par- 
ticle i. In the case of monochromatic radiation, there is an ana- 
lytic solution for photoionisation rate as a function of dept h into 



A4i 



the slab which can be used to determine T| (see Appendix 
For polychromatic spectra we tabulate the solution as a function of 
four variables, {T^ /n n , r far , T, rS car }Q 

The first variable, r far /n H , is the ratio of the amplitude of the 
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Figure 1, An illustration of two density field discretisation methods. The 
top panel shows a multi-valued relationship between points in space and 
resolution elements (SPH). The bottom panel shows a single valued rela- 
tionship between points in space and resolution elements (rectilinear). We 
show a single ray of length Z ray traced from each resolution element. In 
the SPH case, a change in the neutral fraction of any element changes the 
optical depth calculated along all three rays. In the rectilinear case, chang- 
ing the neutral fraction of element 1 (for example) leaves the optical depth 
calculated along the rays from elements 2 and 3 unchanged. 



incident radiation field over the density of the slab, the second, r ar , 
determines the incident spectrum (i.e. how much the user supplied 
spectrum has been hardened before entering the slab), the third, T, 
determines recombination and collisional ionization rates, and the 
fourth, TH Car is the thickness of the slab. In both the analytic and the 
lookup table case, we label this solution Q and note that the values 
of r f ar , r far , and -rg oar are different for each of the 7V ray rays traced 
from particle i. In summary, the total shielded photoionisation rate 



for particle i is computed in the following way. We trace rays 



k = 1, 



y^far 
1 k 



, Ar ra y and calculate, 
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dv 
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Tfar 
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(12) 



(13) 



(14) 



This whole process is illustrated in Fig. [2] The determination of 
r shld has no dependence on the neutral fraction of particle i or 
any of its neighbors and therefore the solution is numerically sta- 
ble. However, it requires extra computational effort to evaluate the 
function Q and introduces errors due to the lack of interaction be- 
tween different near slabs. In the following section we will quantify 
these errors. 



The optical depth 



iVS"(T is a function of frequency and not 



a scalar like the other parameters, however, the full shape of r tar (^) is 
determined by a single evaluation, for example r (^th)- 



Urchin 



a) 



T-*shld 


near 


y-ifar 


far 


y-4hin 


T n R 


near 

t H 




far 
T 





b) 



far 


near 




near 


far 


1 


2 




3 


4 



c) 




Figure 2. An illustration of our reverse ray tracing technique in SPH density 
fields. Panel a) shows a typical situation in Test 2 with a single ray being 
traced from a particle with density n H and temperature T. The far slab on 
the right indicates the optical depth r far due to far particles whic h attenu- 
ate the optically thin radiation field and determine T far (Eq. |l2J . This, in 
turn, determines the radiation incident on the slab used to model the near 
particles. The density and temperature of the near slab are determined by 
the particle being updated (n H ,T) but its thickness rg ear = A r g oar cr th 
is fixed by all near particles. The solution Q is then used to determine the 
photoionisation rate r shld at the bottom of the slab independent of the neu- 
tral fractions of any of the near particles (Eq. [13). The v ariables T,n m ,y, 
and r shld are then used in an analytic solution (Eq. |A9[ to determine the 
updated x for the particle. Panel b) shows a typical situation in Test 3, in 
which two rays are traced from each particle and a r shld is determined for 
each ray. These are then combined to find the total shielded photoionisation 
rate for the particle (Eq. |14) . In this case, radiation that should be incident 
on slab 2 from the direction of slab 3 is not accounted for by the solution Q 
and vice versa. Panel c) shows a 2-D cartoon of the situations in Test 4 in 
which rays are traced in all directions. A component of the errors in these 
tests is due to an extension of the problem in using Q described for panel 
b). 



5 TESTS AND VERIFICATION 

In this section, we present several tests performed in order to vali- 
date URCHIN. We begin with simple test cases with a known ana- 
lytic solution, and end with more realistic tests that involve gaseous 
galactic halos drawn from a cosmological simulation. For those 
tests in which an analytic solution is not available, we compare 
the URCHIN results to those of a straight-forward numerical solver 
which we will call MS to distinguish it from URCHIN. After we 
have verified MS against analytic solutions, we will simply refer 
to MS as the analytic solution. 
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Figure 3. Test 1 : Equilibrium neutral fraction x as a function of depth z in 
the case of plane-parallel radiation incident onto a slab of hydrogen with 
constant and uniform density n H = 1.5 X 10 — 3 cm — 3 and temperature 
T = 10 4 K. The thick black line is the solution when using the spectrum of 
Haardt & Madau 2001 restricted to photons with energy 1 < hu/hu^ < 4 
(iJF )• The pink line is the MS result for the case of monochromatic pho- 
tons with energy hu = 19.2 eV (the grey approximation of I^ M )\ green 
dots are the corresponding analytic result discussed in the Appendix. The 
excellent agreement demonstrates the accuracy of MS. By construction, 
the grey approximation reproduces both the number density of photons, 
n 7 , and the photoionisation rate, T, of I~ in the optically thin limit. The 
dashed yellow and dashed olive curves represent cases in which monochro- 
matic photons with hv = 13.6 eV were incident and the flux was nor- 
malised to reproduce either T or n 7 , respectively. These monochromatic 
approximations to HMO I give inaccurate results either in the optically thin 
limit, and/or in the position of the ionisation front, illustrating the need to 
properly sample the input spectrum. 



5.1 Test 1: Analytic Slab Solution 

This test involves plane-parallel radiation with flux F incident from 
one side onto a slab of hydrogen gas of thickness L s i a b, uniform 
density 7ih, and fixed uniform temperature T s i a b. The surface of 
the slab is coincident with the z — plane and the bulk extends in 
the z > direction. In the case of monochromatic radiation, this 
problem has an analytic solution which we derive in Appendix A4. 
The purpose of this test is to verify MS and to illustrate the de- 
pendence of the solution on the assumed spectrum of the incoming 
radiation. 

The equilibrium neutral fraction as a function of depth, x(z), 
is obtained with MS by dividing the slab into many thin slices per- 
pendicular to the z-axis. Starting from z — and working down- 
wards, we solve for x in one slice at a time using the following al- 
gorithm: 1) determine the HI optical depth, r = f z x n H dz, above 
the current slice; 2) calculate an attenuated photoionisation rate in 
the slice, V — J°° Faexp(— r)dv; 3) determine x in the slice by 
plugging r into an analytic solution (Eq.[A9j. To avoid errors due 
to finite slice width, we choose the number of slices such that each 
has a total hydrogen optical depth th = Njiath below unity. This 
guarantees they will always be optically thin when considered in- 
dividually. The numerical values for the parameters of this test are, 
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isiab = 200 pkpc, n H = 1.5 x 10 cm (500 times the cosmic 
mean n H at redshift 3), and T s i a b = 10 4 K. We compare the ana- 
lytical solution to that obtained by MS in Fig.|5](pink line labelled 
grey versus green symbols): they agree very well. 

An advantage of the reverse ray tracing approach is the high 
fidelity with which the spectrum of the ionising background can 
be sampled. This is non-trivial in forward ray tracing schemes and 
many rely on a simplified spectrum (see discussion in Mirocha et al. 
|2012| >. We can use MS to quantify the accuracy of such approxima- 
tions. The slab parameters were chosen in order to produce a fully 
neutral region on the far side of the slab and a to tal column density 
N m = 10 20 - 3 crrT 2 when the redshift z = 3 JHaardt & Madau| 
2001 (HM01) UV background is incident. In what follows, we will 
refer to the specific intensity of the HM01 UV background at red- 
shift 3 between the frequencies ^th and 4l/th as I^ M . In Figure |3| 
we compare HM01 with three different monochromatic approxi- 
mations: i) 1 Rydberg photons with a flux that produces the same 
(optically thin) density of ionising photons as I^ M , ii) 1 Rydberg 
photons with a flux that produces the same (optically thin) pho- 
toionisation rate as /„ , and Hi) the grey approximation of /Jr 
which reproduces both the number density of photons and pho- 
toionisation rate of /J IM . Photons in the grey approximation have 
an energy of 19.2 eV. 

In all cases the monochromatic solutions transition to fully 
neutral more abruptly than the true HM01 solution. This is due 
to the photons having a fixed photoionisation cross-section. The 
spread in frequencies in the HM01 spectrum smooths the transition 
from highly ionised to fully neutral. In the worst case, this can cause 
an error of several orders of magnitude in the neutral fraction. The 
monochromatic spectrum normalised to the same photo-ionisation 
rate (yellow line) recovers the correct neutral fraction in the opti- 
cally thin region but underestimates the depth of the ionised region 
by almost 100 pkpc. The monochromatic spectrum with the same 
number density of ionising photons (olive line) underestimates the 
neutral fraction in the optically thin region by ~ 0.5 dex, but re- 
covers the location of the ionisation front to better than 10 kpc. The 
grey approximation (pink line) recovers the ionised fraction in the 
optically thin region and the location of the ionisation front, but still 
has a maximum error of ~ 0.75 dex. 

The dependence of x(z) on the spectrum of radiation illus- 
trates the need for accurate multi-frequency treatments. An advan- 
tage of URCHIN is that it can treat arbitrary spectral shapes accu- 
rately with approximately 100 frequency bins. In addition, the full 
attenuated ionising spectrum at a particle position is known when 
neutral fractions are calculated. For completeness we examine sev- 
eral other frequently used approximations to the HM01 spectrum 
in Appendix | A2| and various fitting formula used to approximate 
the hydrogen photoionisation cross section in Appendix |Al| 

5.2 Test 2: Uniform Slab - Radiation Incident from One Side 

In this test we compare URCHIN solutions to those of MS. To 
this end, we create a set of uniform slabs with the same geome- 
try as in Test 1, but vary the volume densities nn such that the 
projected HI column densities through the slabs cover the range 
17.5 < iVm/cm -2 < 20.3 (i.e. the range over which self- 
shielding becomes important). To create SPH realizations of the 
uniform density fields we generate glass-like distributions with 
16 3 , 32 3 , and 64 3 particles (hereafter labeled N16, N32, and N64). 
To model the plane parallel radiation in URCHIN, we trace a single 
ray of length L s i a b from each particle towards the surface of the 
slab. We calculate column densities by projecting all SPH particles 
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Figure 4. Test 2: Plane parallel radiation incident from one side onto 
200 pkpc isothermal slabs of uniform density. The spectrum is that of HMO 1 
at redshift 3 integrated between 1 and 4 Rydbergs. The slabs shown have 
n H = 1.5 x 10- 3 cm- 3 x {1.0, 0.96, 0.92, 0.88, 0.84, 0.80, 0.76, 0.72, 
0.68, 0.64, 0.60 } (from black to red) which, at redshift 3, correspond to 
overdensities of A = {500, 480, 460, 440, 420, 400, 380, 360, 340, 320, 
300}. Top panel: equilibrium neutral fraction as a function of depth x (z) for 
the analytic solution (solid lines) and the URCHIN N16 resolution solution 
(dashed lines). Lines are labelled with the corresponding value of the over 
density (n H = n ™ can x A), and the analytic neutral hydrogen column 
density through the slab, logA r g^ a . Bottom panel: the difference in A'hi 
between the URCHIN and analytic solutions, as a function of the analytic 
Njjx. Shown are the results when the HM01 spectrum is used (solid lines) 
and when the grey approximation is used (dashed lines) , with different 
colours referring to different particle resolutions. In all cases, the URCHIN 
reverse ray tracing solution is within 0. 1 dex of the analytic value, with the 
error depending weakly on particle resolution. 
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onto a plane and measure the mean column on a fine grid of 2048 2 
pixels. Similar projections were used in jAltay et al-ipofl) . 

In the top panel of FigureE] we compare AfS solutions (solid 
lines) to those produced by URCHIN (dashed lines) for slabs with 
different densities (different colours). The URCHIN solutions are 
from the lowest (N16) resolution SPH density fields. URCHIN faith- 
fully follows the dependence of neutral fraction on depth for all 
models, including those where the gas becomes mostly neutral. The 
bottom panel of the figure quantifies the errors in neutral hydrogen 
column densities calculated from the SPH realizations (solid lines). 
Errors are below 0.05 dex at low resolution (N16, blue line), and 
improve with increasing resolution. The dashed lines quantify the 
URCHIN errors in the case of monochromatic radiation, for which 
the analytic solution Q does not require the construction of the in- 
terpolation table discussed in Section [4.1.2| The errors are slightly 
larger here as the ionisation front becomes very steep when the ra- 
diation is monochromatic. However, this test demonstrates that our 
method of splitting the optical depth into a contribution from near 
and far particles works well for single ray applications. 

5.3 Test 3: Uniform Slab - Radiation Incident from Two Sides 

This test is identical to Test 2, except we irradiate the slabs from 
both sides. To obtain the analytic solution, we modify MS as fol- 
lows. First we initialise the neutral fractions of all slab slices to 
their optically thin values, then we loop over the slices calculating 
the optical depth both above and below a given slice. These optical 
depths are used to calculate a photoionisation rate and hence a new 
value for the neutral fraction. We continue iterating over the slices 
until the neutral fraction in each slice has converged to one part in 
ten thousand. The solution is symmetric with respect to the centre 
of the slab at z = 100 pkpc. To model the plane parallel radiation 
in URCHIN, we trace two rays from each particle, one in the +z 
direction and one in the —z direction. 

As in the previous test, we compare the URCHIN (dashed 
lines) and analytic (solid lines) solutions for the neutral fraction 
in the top panel of Fig. [5] The goal of this test is to examine the 
accuracy of the near I far split described in section 4.1.3 when 
multiple rays are being used (diagramed in panel b of Figure [31. 
The algorithm introduces errors in the calculation of r ahld because 
each ray is considered independently. To illustrate this point we 
will consider the process of calculating r shld = 2ir(Tf ld + T sMd ) 
for a particle situated in the middle of a slab in which the equi- 
librium neutral fraction doen not form a neutral core (i.e. any slab 
with A < 460). First the +z ray is traced and a F+ r is calculated 



as in Eq. 1 12 1. This is then used as input to calculate F^_ as in 



Eq.l|13[> which represents the photoionisation rate at the bottom of 
the near slab. The error occurs due to the fact that the near slab 
should also be irradiated from the — z direction as all of the gas 
is highly ionised. This leads to an overestimate of the opacity of 
each of the near slabs and in turn to an overestimate of the neutral 
fraction in the slab. The errors are most severe during the transition 
from optically thin slabs to slabs that form a neutral core. In slabs 
that do form a neutral core the error is absent as at least one ray 
always encounters a high optical depth. 

The resulting errors on the column density through the slab 
are shown in the bottom panel of Fig. B] Different colours refer 
to different numerical resolutions, with solid lines representing the 
HM01 radiation field, and dashed lines the grey approximation. The 
URCHIN optical depth is within 0.1 dex of the analytic result for 
columns JV H i < 10 18 cm" 2 or 7V H i > 10 19 ' 5 cm" 2 for the highest 
resolution slab. These column density limits correspond to the cases 
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Figure S. Test 3: Plane parallel radiation incident from two sides onto 
200 pkpc isothermal slabs of uniform density. The spectrum is that of HMO 1 
at redshift 3 integrated between 1 and 4 Rydbergs. The panels are arranged 
as in Figured Top panel: for slabs with A ^ 340 or A ^ 460, errors in 
the ionisation profiles are similar to Test 2. In the intermediate regime both 
rays are important in determining the photoionisation rate at a particle and 
the assumptions made in the solution Q are not valid (see Figurepl. Bottom 
panel: the difference in Nm between the URCHIN and analytic solutions, 
as a function of the analytic TVjji ■ Because the importance of Q is reduced at 
higher resolutions, these errors scale more strongly with particle resolution 
than those in Test 2. 



where the slab is either mostly ionised everywhere, or develops a 
neutral core. In the intermediate regime, URCHIN overestimates the 
neutral hydrogen column density by up to 0.15 dex at N64 reso- 
lution, and 0.25 dex at N16 resolution. The errors in the case of 
monochromatic radiation (dashed lines) are not significantly dif- 
ferent, showing that the URCHIN error is due to the near/far split, 
rather than the implementation of the look-up table for polychro- 
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matic spectra employed in Q (as opposed to the analytic solution 
used for monochromatic spectra). 



5.4 Test 4: Galactic Halos - Uniform UV Background 

In previous sections we applied URCHIN to very simple slab ge- 
ometries for which we could calculate accurate analytic solutions. 
In this section we concentrate on the more realistic case of galac- 
tic halos. In particular we focus on halos extracted from the OWLS 
suite of cosmological galaxy formation simulations (Scha ye et al.| 
2010). This suite consists of a reference model (REF) and more than 



50 variations around that reference model which explore changes to 
sub-grid physics and other parameters. The REF model, which we 
make use of here, included a model for pressure in the numerically 
unresolved cold interstellar medium, star formation, the timed re- 
lease of 11 chemical elements by type I and type II supernovae 
and AGB stars, radiative cooling due to the same elements in the 
presence of the HM01 ionising background, and energetic feedback 
from supernovae. 



5.4.1 Stacking Halos 

To identify these halos we first calculate a friends-of-friends group 
catalogue using the standard algorithm of |Davis et al.| ( |T985[ > run 
on the dark matter particles, and linking baryonic particles to their 
nearest dark matter particle. We then identify bound sub-structures 
(i.e. halos) within these FoF groups using the subfind algorithm 
( Springel e t al.||206T[ |Dolag et al.|2009) and associate these ha- 
los with galaxies. Our goal is to create objects that have approxi- 
mate spherical symmetry - so we can calculate self- shielding an- 
alytically - yet are representative of typical galaxy density pro- 
files encountered in simulations. To this end, we 'stack' haloes 
of similar mass. First, haloes are grouped into mass bins accord- 
ing to their total gas mass, with the lowest mass bin edge being 
Mgas = 1O 815 M /i _1 (100 gas particles) and the bin width equal 
to 0.3 dex. We then centre all haloes in a given bin on the location 
of the most bound particle of that halo. The mass of each particle 
in a stack is then adjusted such that the sum of all particle masses 
in the stack is at the logarithmic center of the mass bin. Finally, the 
smoothing lengths and densities of all particles in a stack are re- 
calculated. These stacks, while still containing some particle noise, 
are now (nearly) spherically symmetric. In Figure [6] we show im- 
ages of each halo stack out to a radius of 100 pkpc. The color scale 
logarithmically covers the range between Nn = 10 17 cm -2 and 
Nh = 10 23 cm~ 2 . We show this figure mainly to give the reader 
a visual impression of how the halo stacks become less spherically 
symmetric in the higher mass bins due to the smaller number of 
halos in each bin. 



5.4.2 Radial Profiles 

In Fig. (ITl we plot radial density profiles of the stacked haloes. In 
each panel there are four black dashed lines. The lower horizontal 
line indicates the cosmic mean value of n H at z = 3, the upper 
horizontal line indicates the star-formation density threshold used 
in REF, the vertical line indicates the gravitational softening length, 
and the diagonal line is an arbitrarily normalised 1/r 2 line to guide 
the eye. Both the total and neutral hydrogen number density are 
shown with shaded regions. The two quantities are equal at small 
radii causing the shaded regions to have a shape similar to the greek 
letter lambda in each panel. 



5.4.3 Total Hydrogen Radial Profiles 

Total hydrogen density profiles were calculated from the SPH 
stacks and are plotted as the upper blue shaded regions. The width 
of the regions correspond to one sigma variations from the mean 
density at a given radius. This is a measure of the deviations from 
spherical symmetry shown in Figure [6] The shaded regions are 
mass weighted averages but for reference we also plot the volume 
weighted average, i.e. the sum of all particle masses in a radial shell 
divided by the volume of the shell, with a green line in the bot- 
tom left panel. The volume weighted density hugs the lower end 
of the mass weighted density. In an SPH distribution with no par- 
ticle noise the two quantities would be equal; however, any clump- 
ing will increase the mass weighted quantity relative to the volume 
weighted quantity. The offset between the two at small radii indi- 
cates the level of particle noise in the spherically symmetric part 
of the SPH stacks while the differences at large radii are due to 
clumps caused by stacking a finite number of halos. Note that this 
noise was mostly absent in the previous tests due to the use of glass 
like distributions. 

For each mass bin we construct two smooth analytic n H pro- 
files by fitting a polynomial through the one sigma variations de- 
scribed above (upper red shaded region). The differences between 
the two are only visible at the larger radii of the more massive bins. 
We use the point where the volume weighted rt H intersects the cos- 
mic mean n H to define a radius for each analytic profile which is 
why the red shaded regions are truncated at smaller radii than the 
blue. This provides us with a perfectly smooth and spherically sym- 
metric approximation to our SPH stacks. In general, the density at 
a given radius is monotonically increasing with halo mass as is the 
radius of the halos. Each halo has a profile close to 1/r 2 at inter- 
mediate radii but becomes steeper at both smaller and larger radii. 



5.4.4 Neutral Hydrogen Radial Profiles 

To calculate neutral hydrogen profiles in spherically symmetric gas 
with J\fS we did the following. Each halo is divided into N B h shells. 
From each shell, we trace Ne rays which sample the azimuthal an- 
gle between and it radians. The photoionisation rate in each shell 
is calculated by summing the contribution from each ray and is then 
used to calculate a new neutral fraction. We loop over all shells and 
iterate until the neutral fraction in each shell has converged to one 
part in ten thousand. We find that our results are numerically con- 
verged when using N s h = 1600 and Ne — 13. We use MS to 
calculate n HI profiles from the —1 a and +lo~ n H profiles, and 
show the results as the lower red shaded regions. 

We also use URCHIN to calculate the neutral fraction of every 
particle in each stack and construct average n HI profiles in the same 
way as we calculated the n H SPH density profiles. These are shown 
as the lower blue shaded areas. When irradiated from outside, each 
halo develops a characteristic ionisation structure with a neutral 
core, a sharp ionization front, and an optically thin region. In the 
core, n H = n HI . When the density drops to logn H « — 2 the neu- 
tral core transitions to a steep ionisation front. At densities between 
logn H = —4 and Iogn H = —5 the ionisation profiles transition 
into the optically thin regime in which the neutral fraction is pro- 
portional to the density, x — an^ /Y . In all cases, the URCHIN so- 
lution overlaps with the analytic solution. The differences between 
the two are largest in the ionization front but match the analytic 
solution in the optically thin and optically thick regimes. This is 
the spherical corrolary to the situation in Test 3 in which the errors 
were largest for intermediate column density slabs. In a spherical 
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Figure 6. Test 4: Projected total hydrogen density for the nine halo mass bins used for constructing the halo stacks. The color scale logarithmically covers 
the range between Nn = 10 17 cm — 2 and JVh = 10 23 cm — 2 . Each image shows a region 200 pkpc on a side and deep enough to include every particle in 
the stack. These images give a sense of how the higher mass halo stacks are less spherically symmetric than the lower mass stacks, as there are fewer of them 
within the simulation volume to average over. 



geometry, even with a neutral core, some rays will sample the prob- 
lematic columns between 10 18 < A?Hi/cm~ 2 < 10 19,5 cm~ 2 . 



5.4.5 Neutral Hydrogen Column Density Profiles 

In Figure IS] we show the neutral hydrogen column density Nm as 
a function of impact parameter for the halo stacks. To calculate 
JVhi profiles in the analytic case we simply integrate lines of sight 
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Figure 7. Test 4: n H and n HI radial profiles for OWLS halo stacks irradiated with a uniform UV-background. Different panels correspond to different halo 
gas mass ranges as indicated in each panel. The lower horizontal dashed line indicates the cosmic mean value of n H at z = 3, the upper horizontal line 
indicates the star-formation density threshold used in the REF OWLS model, the vertical line indicates the gravitational softening length, and the diagonal line 
shows an arbitrarily normalised 1/r 2 relationship to guide the eye. Both the total (n H ) and neutral (n HI ) hydrogen number density are shown with shaded 
regions. The two quantities are equal at small radii causing the shaded regions to have a shape similar to the greek letter lambda in each panel. The SPH 
density fields on which URCHIN was run are shown in blue while our analytic solutions are shown in red. All SPH profiles (blue) are constructed from mass 
weighted averages at a given radii. The width of the shaded regions indicates one sigma variations from the median, and hence are a measure of deviations 
from spherical symmetry in the stacks. In the lower left panel we also indicate the volume weighted average of n H with a green line. The URCHIN solution 
and the analytic solution always overlap. The two solutions are most different in regions where the n HI profile is steepest. 



through the spherically symmetric shells shown in Figure [7] and 
tabulate the results as a function of impact parameter b. We calcu- 
late these profiles for the —la and +la spherical mass profiles, 
and shade the region between the two in red. To calculate the same 
quantity from the SPH distributions we project all particles in a 
stack onto a plane and measure the column density on a fine grid 
of 2048 2 pixels. We then bin these pixels in impact parameter and 



show the one sigma variation around the mean as a shaded blue re- 
gion. For density profiles of the form n m oc r~ n with n > 1/2 
the dominant contribution to iVm along a line of sight comes from 
the smallest radii. Because of this, we can associate lines of sight 
in Figure [8] with the innermost radius they probe. This gives rise 
to a corresponding region in the Nm(b) plots for each of the re- 
gions described in the previous section (neutral core, ionisation 
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Figure 8. Test 4: Nm impact parameter profiles for OWLS halo stacks irradiated with a uniform UV-background. Different panels correspond to different 
halo gas mass ranges as indicated in each panel. The SPH density fields on which URCHIN was run are shown in blue while our analytic solutions are shown 
in red. The SPH profiles (blue) are constructed by projecting all particles in a halo stack onto a plane. The width of the shaded regions indicates one sigma 
variations from the median value at a given impact parameter. The URCHIN solution and the analytic solution always overlap. 



front, and optically thin outskirts). In all mass bins, the URCHIN 
solution overlaps with the analytic solution giving us confidence 
that URCHIN can be used to accurately model neutral hydrogen ab- 
sorbers in a cosmological context. 



6 DISCUSSION AND CONCLUSIONS 

We have presented and described a new publically available ra- 
diative transfer code called URCHIN. It is optimized to model the 
residual neutral hydrogen in the post-reionization Universe and re- 
lies on reverse ray tracing to avoid problems typically associated 



with modelling the UV background. In particular, our implementa- 
tion allows for the preservation of symmetry in the input radiation 
field, a completely uniform sampling of gas resolution elements, a 
high fidelity sampling of spectral features, and some optimisations 
not possible in forward ray tracing schemes. 

We have validated URCHIN in four sets of tests that range from 
comparison to analytic solutions in simple slab geometries to solv- 
ing for the ionization structure in gaseous galactic halos. The errors 
discussed in §5 have two root causes. One is particle resolution, and 
the other is the multi-valued relationship between points in space 
and resolution elements in SPH. The only solution to the first prob- 
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lem is to run simulations with more particles. The second problem 
prompted us to introduce the near I far split and the analytic slab 
solution Q discussed in §4. This solution solves the multi-valued 
problem but is restrictive in its use of plane parallel geometry and 
its treatment of each ray / near slab pair independently (i.e. with- 
out considering the other pairs). It is possible that an approach that 
parameterises a spherical analytic solution using some average of 
the near slab optical depths and some average of the incoming ra- 
diation from all rays as parameters would be more sucessful. In 
addition, the analytic solution presented in Eq. |A22| considers uni- 
form density slabs although there is a density derivative in the fully 
general solution. Estimating dn n jdr along each ray and using this 
as an extra parameter in Q could also prove useful in improving the 
solution. 

In the current implementation, we have shown that solutions 
generated by URCHIN reproduce analytic solutions in the case 
of density fields approximating gaseous galacitc halos. We view 
URCHIN as a first step in a series of incremental improvements to 
the standard treatment of the UV background. Future steps will in- 
clude helium as in Al tay et aL"|(20 08l, proximity zones from inter- 
nal point sources, an equilibrium heating/cooling model, and non- 
equilibrium corrections for ionisation and heating. 
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the form J(v) = Jo 0/Vth)~ a , f or a = 0, 1, 2. The HM01 spec- 
trum has features due to re-emission, such as the Helium Lyman-a 
emission at 3 Rydberg. The effect of using such approximate spec- 
tra as opposed to the full HM01 spectrum is illustrated in the middle 
panel. We calculated the equilibrium neutral fraction as function of 
depth, Xen(z), for plane-parallel radiation entering a slab with uni- 
form hydrogen density n H — 1.5 x 10 -3 cm -3 (500 times the 
cosmic mean riH at z = 3) and T = 10 4 K. Power-law approxi- 
mations to HM01 work reasonably well, as long as the spectra are 
normalised to give the same number density of ionising photons as 
the orginal HM01 spectrum (dashed lines). Normalising the power- 
law spectra to the same photoisation rate does not work well (solid 
lines). The right panel in Figure [AT] compares the differences in 
neutral columns at z — 0.2 pMpc, between the full HM01 spectrum 
and these powerlaw approximations. The largest error occurs when 
the slab fails to become fully neutral, as in the flat spectrum case 
(a = 0). The monochromatic grey approximation used in Section 
|5.1| (pink bar) overestimates the central HI column by ^ 0.04 dex. 
The bottom two bars show the effect of approximating the pho- 
toionisation cross-section by a simple power law versus using the 
more accurate fit from [Verner et al.| p996), and using the HM01 
spectrum up to 10 Rydberg, as opposed to 4 Rydberg as we have 
done up to now. 



APPENDIX A: ANALYTIC HYDROGEN IONIZATION 
SOLUTIONS 

Al Photoionization Cross Section 

The photoionization cross section for hydrogenic atoms in the 
ground state can be expressed analytically as, 

o-th I v th \ 4 exp {4 - [(4 tan" 1 e)/e] } 



Z 2 



/y th y exp^ 



exp(— 27r/e) 



1/ Z 71 2 
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where a& is the fine structure constant, ao is the Bohr radius, and 
Z is the atomic number of the atom. Two fitting formula are com- 
monly used in the literature. The first is a power law form, usually 
accompanied by a citation to Osterbrock ( 1989 ), 
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This simple form can be useful for analytic treatments, but is not as 
accurate as the fit due to V erner et al.| ( |1996l >, 



a — Co (X — 1) X 
hv 



2 0.5P-5.5 



1 + 



y/x/x a j 



X 



E = 0.4298 eV, a = 5.475 X 10* Mb 
x a = 32.88, P = 2.963. 



(A3) 



We explore the errors these approximations introduce into column 
density calculations in the next section. 

A2 Power-law approximations to the Haardt & Madau 2001 
spectrum 



The left panel in Figure |A1| compares the shape of the Haardt & 
Madau (2001) (HM01) spectrum to power-law approximations of 



A3 Time dependent solution of the neutral fraction 

The rate of change of the hydrogen neutral fraction x = n HI /n H 
is determined by the rate of photoionisation F, collisional ionisa- 
tion 7, and recombination a as well as the number density of free 
electrons n c . 



dx 



= — (T + 7n c )x- + an c (l — x) 



(A4) 



If we decompose the free electron number density into a contribu- 
tion from ionised hydrogen and a contribution yn H from all heavier 
elements, we can write: 



n c = (1 — x + y)n H 
Substituting into the previous equation yields: 



dx 



= -[r + 7 (1 - x + y)n n ]x + 
a(l - x + y)n H (1 - x) 



(A5) 



(A6) 



Grouping terms in powers of x, we can write dx/dt in the form of 
a Riccati equation: 

~ = Rx 2 + Qx + P (A7) 

dt 

R = (7 + a)?i H 

Q = - [F + (7 + 2q)ti h + ( 7 + a)n H y] 

= - [T + qti h + R(l + y)\ 

P = cm H (l + j/). (A8) 

The roots of the quadratic term are 

-Q - (Q 2 - 4PR) 1/2 

X - ~ 2R 

-Q + (Q 2 -4PJ^) 1/2 
X+ 2R 

To determine which of these roots is the physical equilibrium solu- 
tion we consider the case of pure hydrogen (y — 0) in the absence 
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Figure Al. Left panel: The 2: = 3 |Haardt & Madau|{2001} UV-background (HMOl) spectrum compared to three different power-law spectra, J{y) = 
Jo (v/v t h)~ a , for a = 0, 1, 2 shown blue, green and red, respectively. For the solid lines, Jo is chosen such that the power-law spectrum has the same 
photoionisation rate as the HMOl spectrum, for the dashed lines Jo is chosen such that they have the same number density of ionising photons. The spike 
around 3 Rydbergs is due to Helium Lyman-a emission. Middle panel: equilibrium neutral fraction, x Gq (z) as a function of depth z into a slab with uniform 
hydrogen density n H = 1.5 X 10 -3 cm -3 (500 times the cosmic mean «h at 2 = 3) and T = 10 4 K. The thick black line is for the full HMOl spectrum, 
blue, green and red are the corresponding power-law approximations. Such power-law models work relatively well, provided the amplitude Jo of the spectrum 
is chosen such that the spectra have the same number density of ionising photons (dashed lines). Right panel: difference between the total Nm for these 
spectral approximations and that in the HMOl model. For the four bars that go off the plot, we have indicated the horizontal value in the panel. The error in 
Nui remains below 0.05 dex for all cases except those in which the approximating spectra are normalized to have the same optically thin photoionization rate 
in which case the errors can be much larger. Of course the neutral fraction in the optically thin case is incorrect unless spectra are normalised to have the same 
optically thin photoionization rate. For reference, we also show as a brown bar differences in the hydrogen columns for the HMOl spectrum due to using a 
simple power law for the photoionisation cross section, as opposed to using the fit from Verner et al. 1 1996), see Eq.JA3l. Finally the black bar shows the effect 
of truncating the HMOl spectrum at 10 Rydberg, as opposed to 4 Rydberg. 
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The collisional ionisation and recombination rates depend on tem- 
perature T (see for example the fits in Theuns et al. 1998b I which, 
in general, will change as x changes. However in the case of con- 
stant T, n H , y, and T, the coefficients P, Q and R are constants as 
well. We can rewrite the derivative using Vieta's formula as, 

dx 

— — R(x - x cq )(x - x+) (A17) 

and solve for the time-dependent solution by separation of variables 

dx 



X(t) = 

trc =: 



R{x- 

cq + {X0 — 

(x+ 



l )(x-x+) 



dt 



(A 18) 



(a 



F 



cxp 



Cq; (*+ 

3-cq/'' 



xo) + (x — x cq )F 



tr 



(a + 7)n H 



(A 19) 



In a gas composed only of hydrogen and helium, y is bound be- 
tween and (1 — X)/{2X) where X is the hydrogen mass fraction 
of the gas. The solution x(t) is fully determined once values for T, 



n m , y, and T are specified. In addition, all of the dependence on the 
ionisation state of elements other than hydrogen is contained in the 
variable y. For all of the tests performed in this paper we set y — 0. 



A4 Neutral fraction in a plane parallel slab 

In the case of plane parallel monochromatic radiation with a photon 
flux F incident on a semi-infinite slab of hydrogen gas with con- 
stant density and temperature, the ionisation structure x(z) can be 
calculated analytically. We orient our coordinate system such that 
the surface of the slab is coincident with the x — y plane and the 
positive z-axis extends into the slab. The photoionisation rate at z 
is r(z) = Fae~ T ^ z ' — Toe~ T ^ z \ where To is the photoionisa- 
tion rate at z = and a is the photoionisation cross-section for 
the monochromatic radiation. In equilibrium, Eq. l|A6} with y — 
implies: 



r 

ro~ 



Exp 



n„xods 



ToaT 



(A20) 



where we have defined X = a(l — x) — 7(1 — x)x. Taking the 
logarithm of the last two terms and differentiating with respect to z 
gives: 



1 dn H d In X d In Yqx 1 
tjj dz dz dz ) 



(A21) 



Setting the derivative of n H to zero and integrating both sides over 
the interval [0, z] yields the depth Nh<j — zn H a at which the neu- 
tral fraction is x. We can write this inverse solution in terms of the 
equilibrium neutral fraction in the absence of radiation (i.e. in col- 
lisional ionisation equilibrium), a; cc = a/ (a + 7) = x(z — > 00), 
and the value of x at the surface of the slab x(z = 0) = xq. The 
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solution z(x) is then 
1 _ 1 

Xq X 



+ In 



~x(l 


- Xo) 


+ 


-In 


l-x)_ 


.ZoO^ce 


-x)\ 



(A22) 



This allows x(z) to be mapped out. Once x is known, a photoioni- 
sation rate can be determined using Eq. |A20| We use this solution 
to verify our numerical solver (ATS in the text) and then use ATS 
in more general (non-monochromatic) cases to verify URCHIN. In 
addition, Eq. |A22j forms the basis of the solution Q. 



